module input_output
  use parameters
  implicit none
  
contains
  
  ! 读取输入参数
  subroutine green_input()
    character(len=80) :: list
    real(kind=8) :: mus, shift, x0, y0, xs, ys
    integer :: i, j
    
    ! 打开输入文件
    open(10, file='grt.dat', form='formatted', status='old')
    
    ! 读取接收器参数
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') list  ! 跳过标题行
    read(10, *) x0    ! 接收器x坐标
    read(10, *) y0    ! 接收器y坐标
    read(10, *) source%z0  ! 接收器深度
    
    ! 读取源参数
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') list  ! 跳过标题行
    read(10, *) source%M0  ! 标量矩
    read(10, *) source%mino, source%maxo  ! 频率范围
    read(10, *) source%no  ! 频率点数
    read(10, *) source%minv, source%maxv  ! 速度范围
    read(10, *) source%nv  ! 速度点数
    read(10, *) source%rake, source%dip, source%str  ! 断层参数
    read(10, *) xs  ! 源x坐标
    read(10, *) ys  ! 源y坐标
    read(10, *) source%zs  ! 源深度
    
    ! 计算源-接收器距离和方位角
    source%r0 = sqrt((x0-xs)**2 + (y0-ys)**2)
    if (y0 >= ys) then
      source%fai0 = acos((x0-xs)/source%r0)
    else
      source%fai0 = pi2 - acos((x0-xs)/source%r0)
    end if
    
    ! 读取计算参数
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') list  ! 跳过标题行
    read(10, *) comp%m  ! 采样参数
    read(10, *) comp%Twin  ! 时间窗口
    
    ! 读取源类型
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') SourceType  ! 源类型
    
    ! 读取源参数
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') Type_STF  ! STF类型
    read(10, '(A)') list  ! 跳过标题行
    read(10, *) source%tou  ! 上升时间
    read(10, *) shift  ! 偏移
    read(10, *) source%fc  ! 中心频率
    
    ! 读取窗口参数
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') WinSwitch  ! 窗口开关
    read(10, '(A)') WinType  ! 窗口类型
    read(10, *) f1, f2, f3, f4  ! 频率参数
    
    ! 读取输出参数
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') Outname  ! 输出文件名
    
    ! 读取层参数
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') list  ! 跳过标题行
    read(10, '(A)') list  ! 跳过标题行
    read(10, *) nly  ! 层数
    read(10, '(A)') list  ! 跳过标题行
    
    ! 初始化媒质参数
    call init_media_parameters()
    
    ! 读取各层参数
    do i = 1, nly
      read(10, *) j, media%z(i-1), media%rho(i), media%vs0(i), media%vp0(i), media%qs(i), media%qp(i)
    end do
    
    close(10)
    
    ! 角度单位转换
    source%dip = source%dip * pi / 180.0d0
    source%rake = source%rake * pi / 180.0d0
    source%str = source%str * pi / 180.0d0
    
    ! 标量矩转换
    if (source%M0 <= 12) then
      ! 当M0小于12时，M0应为地震震级而非矩
      ! 根据Kanamori方程转换为标量矩(dyne cm)
      source%M0 = 10**(1.5*source%M0 - 5.9)
    else
      source%M0 = source%M0 * 1d-22
    end if
    
    ! 确定源和接收器所在的层
    ls = layernumb(source%zs, nly, media%z)
    lo = layernumb(source%z0, nly, media%z)
    print *, 'ls=', ls, '     lo=', lo
    print *, source%zs,source%z0
    print *, media%z
    mus = (media%vs0(ls)**2) * media%rho(ls)
    
    ! 检查输入参数
    if (nly < ls .or. ls < 1) then
      print*, 'Error: info@input-checking_1'
      stop
    end if
    
    ! 如果源在最后一层，添加一个虚拟界面
    if (nly == ls) then
      if (source%zs > media%z(nly-1)) then
        print *, 'Add one fictious interface'
        media%z(nly) = source%zs + 5.0
        media%vs0(nly+1) = media%vs0(nly)
        media%vp0(nly+1) = media%vp0(nly)
        media%qs(nly+1) = media%qs(nly)
        media%qp(nly+1) = media%qp(nly)
        media%rho(nly+1) = media%rho(nly)
        nly = nly + 1
        print*, 'nly=', nly, '     z(nly-1)=', media%z(nly-1)
      else
        print*, 'Error: info@input-checking_2'
        stop
      end if
    end if
    
    if (nly > NM) then
      print*, 'nly>NM !!!'
      print*, 'Current NM = ', NM, '  Current layer number=', nly
      stop
    end if
  end subroutine green_input
  
  ! 计算层号的函数
  ! function layernumb(z0, n, z) result(layer)
  !   implicit none
  !   integer, intent(in) :: n
  !   real(kind=8), intent(in) :: z0
  !   real(kind=8), intent(in) :: z(0:n)
  !   integer :: layer
  !   integer :: i
  
  !   layer = 1
  !   do i = 1, n
  !     print *,i, 'z0=', z0, '     z(i)=', z(i), 'layer=', layer
  !     if (z0 < z(i)) then
  !       layer = i
  !       exit
  !     end if
  !   end do
  ! end function layernumb
  function layernumb(z0,n,z)
    !c       ........
      real*8    z0,z(0:n)
      integer layernumb, n, i
      do i = n-1,0,-1
          if (z0.gt.z(i)) then
            layernumb = i+1
            return
          end if
      end do
      layernumb=1
    !c       print*,'Error: info@layernumb'  
    
    return
    end

end module input_output